* *************************************
* COX PROPORTIONAL HAZARD MODELS
* *************************************
stset year, id(citynum) failure(ownership==1) scale(1) origin(ownership==0) enter(time 1996)
drop if mi(ownership)
*****Core Model*****
stcox title_apply newtypebus brt ridershipyrwan newbuspercent, vce(robust)
****Main model****
stcox title_apply newtypebus brt ridershipyrwan newbuspercent busper1000 gdppercapwan fishealth mayortenure sectenure, vce(robust) nohr
****Main model robustness checks****
stcox title_apply newtypebus brt ridershipbypop newbuspercent busper1000 gdppercapwan fishealth mayortenure sectenure, vce(robust) nohr
stcox title_apply newtypebus brt ridershipyrwan newbuspercent busper1000 gdppercapwan fishealth mayortenure sectenure, vce(robust) nohr
stcox title_apply newtypebus brt ridershipyrwan newbuspercent busper1000 gdppercapwan fishealth mayortenure sectenure, vce(cluster citynum) nohr

*****check for proportionality
estat phtest, detail

* ************ 
* plot residuals
* *************
. predict cs, csnell
. stset cs, failure(ownership)
. sts generate H = na
. line H cs cs, sort xlab(0 1 to 4) ylab(0 1 to 4)
ltable year ownership, survival tvid(citynum)

* ******************
* Dosage frequency effects*****
* **************************
stcox i.history3yr busper1000 newbuspercent ridershipyrwan gdppercapwan fishealth mayortenure sectenure, vce (robust)
stcox i.history5yr busper1000 newbuspercent ridershipyrwan gdppercapwan fishealth mayortenure sectenure, vce (robust)

* ******************
*Test the wellness of fit of Guangzhou case using any of the Cox models****
* ******************
predict dev, deviance
gen dev_abs = abs(dev) if _d == 1
keep if _d == 1
summarize dev_abs if citynum == 196 & year == 2007
scalar focal_dev = r(mean)
twoway    (kdensity dev_abs, lwidth(medthick)) , xline(`=focal_dev', lcolor(red) lwidth(thick))       

**********Weibull survival model*********
streg title_apply newtypebus brt busper1000 newbuspercent ridershipyrwan gdppercapwan fishealth mayortenure sectenure, dist(weibull) vce(robust)

*********** Gompertz survival model***********
streg title_apply newtypebus brt busper1000 newbuspercent ridershipyrwan gdppercapwan fishealth mayortenure sectenure, dist(gompertz) vce(robust)

*********** Discrete-Time Hazard Model (Complementary Log-Log) **********
gen fail = (ownership == 1)
cloglog fail i.year title_apply newtypebus brt  busper1000 newbuspercent ridershipyrwan gdppercapwan fishealth mayortenure sectenure, vce(cluster citynum)


